## Agroforestry Paper Analysis
## Author : l.orero@cgiar.org
## This is an script to analyze and create tables and graphs for paper
## Last update: 27th May 2021

# Create working datasets 
###########################################################################
# import tree dataset
library(readxl)
nyando_1 <- read_excel("2019_Nyando_Data.xlsx", sheet = "set")
View(nyando_1)

# Clean column names
library(janitor)
nyando_1 <- clean_names(nyando_1)

# Convert block and region to categorical variables

nyando_1$block = as.factor(nyando_1$block)
nyando_1$type_r = as.factor(nyando_1$type_r)

# Create subsets
########################################################################

# Create subset for Total Lower Project
nyando_lp <- subset(nyando_1, type_r == "Lower Project")

# Create subset for Total Lower  Control
nyando_lc <- subset(nyando_1, type_r == "Lower Control")

# Create subset for Total Middle  Project
nyando_mp <- subset(nyando_1, type_r == "Middle Project")

# Create subset for Total Middle  Control
nyando_mc <- subset(nyando_1, type_r == "Middle Control")

# Create a population DB for all data
##########################################################################
library(tidyverse)

# create copy of file
nyando_count <- nyando_1
# add new column called count
nyando_count$count <- 1
# select specific columns
nyando_count_1 <- nyando_count %>% 
    select(interviewee_name, block, type_r, area_hectare, count) 

# aggregate these    
nyando_count_2 <- aggregate(. ~ interviewee_name + block + type_r,
                            nyando_count_1, sum)

# Create 5 DBH Classes of five DBH classes of [<10, 10-20, 20-30,
# 30-40, and >40 cm]

#########################################################################

## Tree density of less than 10 cm (<10cm)
nyando_dbh_10 <- filter(nyando_1, dbh < 10)
nyando_dbh_10$count <- 1
## filter required columns
nyando_dbh <- nyando_dbh_10 %>% 
    select(interviewee_name, block, type_r, area_hectare, count) 
## aggregate data for 10cm
aggdata_1 <- aggregate(. ~ interviewee_name + block + type_r,
                       nyando_dbh, sum)
## look up land size from previous aggregated data
aggdata_1$area_hectare <- nyando_count_2$area_hectare[match(
    aggdata_1$interviewee_name, nyando_count_2$interviewee_name)]
## add tree density column
aggdata_1 <- aggdata_1 %>% 
    mutate(tree_d = count/area_hectare)

# Check equality of variance
res <- var.test(tree_d ~ block, data = aggdata_1)
res

##Note: p<0.01 there is significant  difference between the two variances

oneway.test(tree_d ~ block, data = aggdata_1)  #ANOVA Unequal var

nyando_count_2 %>% filter(area_hectare == 0)

## One-way analysis of means (not assuming equal variances) TMP vs TMC
oneway.test(tree_d ~ type_r, data = subset(aggdata_1, 
    block == "Middle"))

## One-way analysis of means (not assuming equal variances)

oneway.test(tree_d ~ type_r, data = subset(aggdata1, 
    Block == "Lower"))
## One-way analysis of means (not
## assuming equal variances)

# 10-20cm
nyando_d2 <- filter(nyando_1, DBH >= 
    10 & DBH <= 20)
nyando_d2$count <- 1
nyando_d21 <- nyando_d2[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]
aggdata2 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d21, sum)
aggdata2$ha <- nyando_count2$hectare[match(aggdata2$IntervieweeName, 
    nyando_count2$IntervieweeName)]
aggdata2$tree_d <- aggdata2$count/aggdata2$ha

## Lower vs Middle ANOVA
aov.2 <- aov(tree_d ~ Block, data = aggdata2)
summary(aov.2)

oneway.test(tree_d ~ Block, data = aggdata2)  #ANOVA Unequal var

oneway.test(tree_d ~ type_r, data = subset(aggdata2, 
    Block == "Middle"))

oneway.test(tree_d ~ type_r, data = subset(aggdata2, 
    Block == "Lower"))

# 20-30CM
nyando_d3 <- filter(nyando_1, DBH > 20 & 
    DBH <= 30)
nyando_d3$count <- 1
nyando_d31 <- nyando_d3[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdata3 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d31, sum)
aggdata3$ha <- nyando_count2$hectare[match(aggdata3$IntervieweeName, 
    nyando_count2$IntervieweeName)]

aggdata3$tree_d <- aggdata3$count/aggdata3$ha


## Lower vs Middle ANOVA
aov.3 <- aov(tree_d ~ Block, data = aggdata3)
summary(aov.3)

oneway.test(tree_d ~ Block, data = aggdata3)  #ANOVA Unequal var

## Control vs Project ANOVA
aov.3b <- aov(tree_d ~ type_r, data = aggdata3)
summary(aov.3b)

# TukeyHSD
TukeyHSD(aov.3b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata3, 
    Block == "Middle"))

oneway.test(tree_d ~ type_r, data = subset(aggdata3, 
    Block == "Lower"))

# 30-40cm
nyando_d4 <- filter(nyando_1, DBH > 30 & 
    DBH <= 40)
nyando_d4$count <- 1
nyando_d41 <- nyando_d4[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdata4 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d41, sum)
aggdata4$ha <- nyando_count2$hectare[match(aggdata4$IntervieweeName, 
    nyando_count2$IntervieweeName)]

aggdata4$tree_d <- aggdata4$count/aggdata4$ha

## Lower vs Middle ANOVA
aov.4 <- aov(tree_d ~ Block, data = aggdata4)
summary(aov.4)

oneway.test(tree_d ~ Block, data = aggdata4)  #ANOVA Unequal var
## 

## Control vs Project ANOVA
aov.4b <- aov(tree_d ~ type_r, data = aggdata4)
summary(aov.4b)

# TukeyHSD
TukeyHSD(aov.4b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata4, 
    Block == "Middle"))

oneway.test(tree_d ~ type_r, data = subset(aggdata4, 
    Block == "Lower"))
## 

# >40cm
nyando_d5 <- filter(nyando_1, DBH > 40)
nyando_d5$count <- 1
nyando_d51 <- nyando_d5[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdata5 <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_d51, sum)
aggdata5$ha <- nyando_count2$hectare[match(aggdata5$IntervieweeName, 
    nyando_count2$IntervieweeName)]

aggdata5$tree_d <- aggdata5$count/aggdata5$ha

## Lower vs Middle ANOVA
aov.5 <- aov(tree_d ~ Block, data = aggdata5)
summary(aov.5)

oneway.test(tree_d ~ Block, data = aggdata5)  #ANOVA Unequal var
## 

## Control vs Project ANOVA
aov.5b <- aov(tree_d ~ type_r, data = aggdata5)
summary(aov.5b)

# TukeyHSD
TukeyHSD(aov.5b)

## Control vs Project ANOVA (Relaxing
## the homogeneity of variance
## assumption)
oneway.test(tree_d ~ type_r, data = subset(aggdata5, 
    Block == "Middle"))
## 

# One way ANOVA DBH1 Lower vs Middle
db1 <- aov(DBH ~ type_r, data = nyando_1_d1)
summary(db1)

nyando_c <- nyando_1
nyando_c$count <- 1
nyando_c1 <- nyando_c[, c("IntervieweeName", 
    "Block", "type_r", "hectare", "count")]

aggdatac <- aggregate(. ~ IntervieweeName + 
    Block + type_r, nyando_c1, sum)
aggdatac$ha <- nyando_count2$hectare[match(aggdatac$IntervieweeName, 
    nyando_count2$IntervieweeName)]
aggdatac$tree_d <- aggdatac$count/aggdatac$ha

## Lower vs Middle ANOVA
aov.c <- aov(tree_d ~ Block, data = aggdatac)
summary(aov.c)

oneway.test(tree_d ~ Block, data = aggdatac)  #ANOVA Unequal var


## Control vs Project ANOVA
aov.cb <- aov(tree_d ~ type_r, data = aggdatac)
summary(aov.cb)

# TukeyHSD
TukeyHSD(aov.cb)

# Levene's Test for Homogeneity of
# Variance (center = median)
oneway.test(tree_d ~ type_r, data = aggdatac)

# AGB
# Calculations############################################
# AGB Calculations Combined sample
# of all DBH
nyando_gb <- nyando_1
nyando_gb$num <- rep(1:length(nyando_1$IntervieweeName))
nyando_agb <- nyando_gb[, c("Block", 
    "type_r", "Agroforestry Practice", 
    "AGB (Mg)", "num")]
aggdata_gb <- aggregate(. ~ num + Block + 
    type_r + `Agroforestry Practice`, 
    nyando_agb, sum)
aggdata_gb$agb = aggdata_gb$`AGB (Mg)`
aggdata_gb$agfp = aggdata_gb$`Agroforestry Practice`
aggdata_gb$Block = as.factor(aggdata_gb$Block)
aggdata_gb$type_r = as.factor(aggdata_gb$type_r)

## Lower vs Middle ANOVA for AGB
oneway.test(agb ~ Block, data = aggdata_gb)  #ANOVA Unequal var

## Lower Project AGB
aggdata_gb1 <- filter(aggdata_gb, type_r == 
    "Lower Project")

# Total AGB under Boundary Planting
sum(filter(aggdata_gb1, agfp == "Boundary planting")$agb)


# Total AGB under Riparian buffers
sum(filter(aggdata_gb1, agfp == "Riparian buffers")$agb)
## [1] 0.09491625

## Lower Control AGB
aggdata_gb2 <- filter(aggdata_gb, type_r == 
    "Lower Control")

## Middle Project AGB
aggdata_gb3 <- filter(aggdata_gb, type_r == 
    "Middle Project")

# Middle Control AGB
aggdata_gb4 <- filter(aggdata_gb, type_r == 
    "Middle Control")
